In-Class Exercise 9

In this exercise, we will import geospatial data, create buffers for accessibility analysis, use spatial joins to assess facility proximity, visualize spatial data on maps, and calculate Hansen’s Accessibility metrics.

Author

Teng Kok Wai (Walter)

Published

October 28, 2024

Modified

October 31, 2024

1 Exercise Reference

2 Overview

In this exercise, we will import geospatial data, create buffers for accessibility analysis, use spatial joins to assess facility proximity, visualize spatial data on maps, and calculate Hansen’s Accessibility metrics.

3 Import the R Packages

The following R packages will be used in this exercise:

Package Purpose Use Case in Exercise
tidyverse Data manipulation and visualization in R. Data wrangling, cleaning, and visualization.
sf Handling and visualizing geospatial data. Importing and managing geospatial data, creating buffers.
SpatialAcc Measuring spatial accessibility metrics. Computing Hansen’s Accessibility.
tmap Thematic map visualization. Creating maps for spatial data visualization.
ggstatsplot Enhanced statistical data visualization. Visualizing accessibility differences across regions.

To install and load these packages, use the following code:

pacman::p_load(SpatialAcc, sf, tidyverse, tmap, ggstatsplot)

4 Data Wrangling

We will download and use the CHAS Clinics and Eldercare Services data sets from data.gov.sg portal.

To read the ELDERCARE shapefile data:

eldercare <- st_read(dsn = "data/raw_data",
                     layer = "ELDERCARE") %>%
  st_transform(crs = 3414)
Reading layer `ELDERCARE' from data source
  `/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/raw_data'
  using driver `ESRI Shapefile'
Simple feature collection with 133 features and 18 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21

To read the CHAS clinics kml file:

CHAS <- st_read("data/raw_data/CHASClinics.kml") %>%
  st_transform(crs = 3414)
Reading layer `MOH_CHAS_CLINICS' from data source
  `/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/raw_data/CHASClinics.kml'
  using driver `KML'
Simple feature collection with 1193 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

We will also use the datasets from Hands-on Exercise 9

mpsz <- st_read(dsn = "data/geospatial",
                layer = "MP14_SUBZONE_NO_SEA_PL") %>%
  st_transform(crs = 3414)
Reading layer `MP14_SUBZONE_NO_SEA_PL' from data source
  `/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
hexagons <- st_read(dsn = "data/geospatial",
                   layer = "hexagons") %>%
  st_transform(crs = 3414)
Reading layer `hexagons' from data source
  `/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
  using driver `ESRI Shapefile'
Simple feature collection with 3125 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 21506.33 xmax: 50010.26 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
eldercare <- st_read(dsn = "data/geospatial",
                     layer = "ELDERCARE") %>%
  st_transform(csr = 3414)
Reading layer `ELDERCARE' from data source
  `/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
  using driver `ESRI Shapefile'
Simple feature collection with 120 features and 19 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21 / Singapore TM
ODMatrix <- read_csv("data/aspatial/OD_Matrix.csv", 
                     skip = 0)

4.1 Data Cleaning and Updating Attributes

4.1.1 Supply:

Set to 100 first to calibrate the model.

eldercare <- eldercare %>%
  select(fid, ADDRESSPOS) %>%
  mutate(capacity = 100)

4.1.2 Demand: 

hexagons <- hexagons %>%
  select(fid) %>%
  mutate(demand = 100)

4.1.3 OD Matrix 

distmat <- ODMatrix %>%
  select(origin_id, destination_id, total_cost) %>%
  spread(destination_id, total_cost)%>%
  select(c(-c('origin_id')))

Convert the distance matrix to km units:

distmat_km <- as.matrix(distmat/1000)

5 Count Number of CHAS Clinics within a Distance Around Each Eldercare Centre

To count the number of points within a distance, we will do the following steps:

  1. create a buffer of 1 km around each eldercare centre
  2. visualize the data
  3. count points

5.1 Create Buffer

Note the singapore crs is in metres. To create a 1km buffer, dist should be 1000. 

buffer_1km <- st_buffer(eldercare, 
                        dist = 1000)

5.2 Visualize Plots 

tmap_mode("view")
tm_shape(buffer_1km) +
  tm_polygons() +
tm_shape(CHAS) +
  tm_dots()

5.3 Count Points

To count number of CHAS Clinics within 1 KM around each Eldercare Centre 

buffer_1km$pts_count <- lengths(
  st_intersects(buffer_1km, CHAS))
buffer_1km$pts_count
  [1] 17 31 26 13 26 28 13 24 11 14 17 32  6 20 21 32 15 19 23 19 11 17 20 22 41
 [26] 17 18 39 24 25 15 16 28 28 14 30 23  6 21 15  9 24 24 25 30  8 32 22 17 28
 [51] 33 16 16 16  8 23 16 12 17 19  8 21 14 20 12 22 11 24 18 17 13 16 13 18 19
 [76] 12 19 23 23 21 15 10 15 23  9 22 26 28 22 16 14 31 25 30  9 14 14 14 10 22
[101] 34 19 16 11  9 17 24 12 31 28 17 16 26 31 24 32 22 29 18 30

6 Computing Hansen’s Accessibility

To compute Hansen’s Accessibility:

acc_Hansen <- data.frame(ac(hexagons$demand,
                            eldercare$capacity,
                            distmat_km, 
                            #d0 = 50,
                            power = 2, 
                            family = "Hansen"))

Then we tidy the output:

colnames(acc_Hansen) <- "accHansen"

acc_Hansen <- as_tibble(acc_Hansen)

hexagon_Hansen <- bind_cols(hexagons, acc_Hansen)

6.1 Visualising Hansen’s Accessibility 

mapex <- st_bbox(hexagons)

tmap_mode("plot")
tm_shape(hexagon_Hansen,
         bbox = mapex) + 
  tm_fill(col = "accHansen",
          n = 10,
          style = "quantile",
          border.col = "black",
          border.lwd = 1) +
tm_shape(eldercare) +
  tm_symbols(size = 0.1) +
  tm_layout(main.title = "Accessibility to eldercare: Hansen method",
            main.title.position = "center",
            main.title.size = 2,
            legend.outside = FALSE,
            legend.height = 0.45, 
            legend.width = 3.0,
            legend.format = list(digits = 6),
            legend.position = c("right", "top"),
            frame = TRUE) +
  tm_compass(type="8star", size = 2) +
  tm_scale_bar(width = 0.15) +
  tm_grid(lwd = 0.1, alpha = 0.5)

6.2 Visualize Statistical Graphic 

hexagon_Hansen <- st_join(hexagon_Hansen, mpsz, 
                          join = st_intersects)
ggbetweenstats(
  data = hexagon_Hansen,
  x = REGION_N,
  y = accHansen,
  type = "p")